Tiltle: Closed-form and numerically-stable solutions to problems related to 
the optimal two-impulse transfer between specified terminal states of Keplerian 
orbits 

The first part of the paper presents some closed-form solutions to the optimal 
two-impulse transfer between fixed position and velocity vectors on Keplerian 
orbits when some constraints are imposed on the magnitude of the initial and 
final impulses. Additionally, a numerically-stable gradient-free algorithm with 
guaranteed convergence is presented for the minimum delta-v two-impulse trans- 
fer. In the second part of the paper, cooperative bargaining theory is used to 
solve some two-impulse transfer problems when the initial and final impulses 
are carried by different vehicles or when the goal is to minimize the delta-v and 
the time-of-flight at the same time. 
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1 Introduction 

The optimal two-impulse transfer between fixed position and velocity vectors on 
Keplerian orbits problem has been studied intensively in the past. Numerical 
and closed-form solutions to various formulations of the the minimum delta-v 
problem have been presented in the literature (see for example [1] ,[?] , [6], [7] or 
[8]). In this paper the optimization problem is formulated in a different way. 
If we considered the case where the 2 impulse maneuvers are implemented by 
two different vehicles (e.g. Earth-departure stage (EDS) and a crew exploration 
vehicle (CEV)) and considering that the unused A Veds cannot be re-utilized, 
minimizing J = Aveds + Avcev might not lead to the optimal performance. If 
we take into account that the Earth-departure stage can provide a maximum Av 
(in interplanetary trajectories this would be a maximum Voo), we can formulate 
the problem in a different way. Given the initial and final states: ro,Vo and 
r f,Vf, calculate two impulse maneuvers such that: 


J = min Avcev 

subject to A veds < Av ma x 

In this paper, we show that a closed-form solution of the above problem 
can be obtained by solving a A th order polynomial. A similar solution can 
be also obtained for the cases when a constraint in Avcev is added to the 
previous problem or when the goal is to minimize the first maneuver subject to 
a constraint in the magnitude of the second maneuver. 

When the optimization index is 


J = Ai>i + Av2 


( 1 ) 
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no closed-form solution can be obtained according to ([?] and [4]). In this 
paper we present a numerically-stable gradient-free algorithm with guaranteed 
convergence. The solution consists of solving two 4 th order polynomials and 
then using the bisection method developed in[5]. 

The final sections of the paper are devoted to solve the previous problems 
using Cooperative Bargaining theory (see [3]). If the two impulsive maneuvers 
are implemented by two different vehicles, we can formulate the problem as 
two agents trying to obtain their maximum benefit (minimizing its Aw) while 
cooperating with a second agent. In Cooperative Bargaining theory Eq. 1 is 
known as the Utilitarian solution, (see [3]). In this paper, we examine some 
properties of the Utilitarian solution (the most common one used in trajectory 
optimization problems) and also the properties of some other solutions known 
in the Cooperative Bargaining literature that could be used to solve the two- 
impulse transfer problem. Once the properties of the Cooperative bargaining 
solutions have been explained, an example of the application of this theory to 
the problem of minimizing the time of flight of the transfer while minimizing 
the total Aw is presented. 


2 Problem setup 

The problem we are trying to solve consists on calculating a two-impulse trans- 
fer between two known states: r^v.^o, with i = 1,2 defining the departure 
and arrival states respectively. In order to solve the optimization problem, the 
following the formulation will be used (see also Battin): 


where: 


vi = v c + v p i v 2 = v c - v p2 


v r = v r xi r 


r i 

Ti 



c = \\r 2 -ri|| 


Cy/jL _ — cos 6) 

rqrqsin# p sin# 


9 £ [0, 7r] angle between iq and r 2 , x = ±i/p, p : parameter of transfer conic 

The sign of x will determine the type of transfer: short (< ir) or long ( n). 
No multi-rev solutions will be considered in this paper. 

It is well known that the velocity locus for will correspond to two hy- 
perbolas (see Sun, Battin) for the short and long transfers (see Figure 1). A 
complete description of the parameters of the hyperbola can be found in Sun 
“on optimum..”. The parameters used in this paper are described in Table 1. 
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Figure 1: Geometry of the hyperbolic velocity locus for terminal 1. Since the 
initial velocity vpo is inside the evolute, four extrema are found. Three of them 
when x > 0. A m j n i(^), , A m j n 2 (^) and one for x < 0; 


Table 1: Geometrical Elements of the hyperbolic locus of velocity 


Asymptotes 

iiri? i c 


Included between asymptotes 

<t> 

Semi-transversal axis 

A h = 2y/v c v p cos \ = ^2 tan § cot \ 

Semi-conjugate axis 

B h = 2 y/v c v p sin § 
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Table 2: Characterization of the extrema of the solutions to the minimum A Wj 
transfer problem. 


1 local minimum for x > 0 

1 local minimum for x < 0 

2 local minima, 1 local maximum for x > 0 

1 local minimum for x < 0 

1 local minimum for x > 0 

2 local minima, 1 local maximum for x < 0 


3 Computing the optimal Avi$ 

A tangent vector to the velocity locus is defined by (see Battin): 
v t i := v c - v p i v t2 := v c + v p2 

Given the departure or arrival velocities v.j ; o we can calculate the parameter 
x that minimizes the departure or arrival A Vi by simply: 

(v c ± v pi ) T Avi = 0 (2) 

Alternatively, we can obtain the same result by 

min Avi = min ||vj - v i>0 || 


and solving for x in 


d ||v» ~ Vi,o|| = Q 

dx 

The solution to eq. 2 will require the solution of the following quartic poly- 
nomials on x for terminals 1 and 2 respectively: 

V C X 4 - V c (i c • Vi o) x 3 + v p (i r i • vi, 0 ) x — w p = 0 (3) 

v c x 4 - v c (i c • v 2j0 ) x 3 - v p (i r2 • v 2j0 ) x - w p = 0 (4) 

Given the signs of the coefficients of the polynomials in eqs. 3 and 4, the 
number of solutions that we can find for x are described in table **. A similar 
approach is described in Sun and other. The solutions provide by both equations 
will create an envelope(s) of trajectories with very interesting properties. For a 
complete description of the trajectory envelopes see Sun. 

We can now compute the global optimal by simply computing the Av associ- 
ated to the real solutions of the quartic polynomials. Let ar,-be the real solutions 
to Eq 3 or 4. Let Av(xj) be the associated Aw, the global optimal solution 

Aw* = min Aw (aw,) (5) 

x j 
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3.1 Existence of multiple extrema 

It is worth studying the case when multiple extrema exist. In Figure ** we can 
see an example of such a case. In this example, the initial velocity vector v 10 is 
inside the evolute associated with the velocity locus of vi. Therefore, there are 
three normals to the hyperbola that intersect at Vyo, that is why three extrema 
appear in the solution eq. 3. The normal in the middle will correspond to a local 
maxima while the two other extreme will correspond with two local minima. 

Given a specific problem, we can calculate a necessary condition for the 
magnitude of the velocity vectors Vi o such that v, j0 is inside the hyperbola 
evolute. For that purpose we will compute the magnitude of the velocity at the 
cusp of the evolute 

_ A 2 h + Bl _ 2 sjv c v p _ A h 

Vcusp — A — d>i ~ 2 d>i 

Ah cos ^ COS 2 ^ 

Ah is also the initial velocity of the minimum-energy transfer orbit v m i, and 
it always corresponds with an elliptical orbit: 



2 tan - cot — 
2 2 


V2 


2 tan 


0 

2 


cot 



Therefore, even initial (or final) elliptical parking orbits could be such that 
Vj Q is inside the hyperbola evolute. 

Finally, a necessary condition for v;.o to be inside the evolute is: 


Vi , o > 


cos 2 ^ 


( 6 ) 


4 Computing the optimal Av\^ when constraints 
in x are present 

In the previous section we have seen that in order to compute the optimal 
departure or arrival A Vi it is necessary to solve a a quartic polynomial in x. 

We have also seen that given the geometry of the problem several local optima 
can appear. If no constraints are present, we can use Eq. 5 to compute the 
global solution. But if constraints on p (or x ) are present a different approach 
is necessary. 

Let x £ C x = [x c i,x c 2 ] be the constraint set for x, with 0 < x c i < x C 2 . An 
algorithm to compute the global optimal Aiq can be found (a similar one Av 2 
can be obtained). 

1. Solve Eq. 3. 

2a. For the real positive solutions of 3. If only one extrema is found, x m i n , get 
Aw* from Table 3. If three extrema are found x m i n \ < x max < x m i n 2 , 
s x = [x m i n i , Xmin?] , get Av* from Table 4 
2b. Repeat 2a for the real negative solutions of 3 
3. The global Aw* will the minimum between the steps 2a and 2b. 
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Table 3: Algorithm for the one-extrema case 

mini C x [x cU X c2 ] 


min ^ C x 

* 

<1 

F 

if C x < X min Aw* = Av(x c2 ) 
if C x > Xmin Aw* = Aw(a: c i) 

T 

Aw* = Av(x min ) 


Table 4: Algorithm for the three-extrema case 


%minl ^ %max ^ %min2i [%minl 5 %min2\i C x 


%minl ^ C x 

%max ^ C x 

-Emin2 ^ C x 

Aw* 

F 

F 

F 

if C x < S x — > Aw* = Aw(a: C 2) 
if C x > S x — > Aw* = Aw(cc c i) 

F 

F 

T 

Aw* = Aw (a : min2 ) 

F 

T 

F 

Aw* = min [Aw(x c i), Aw(a; C 2)] 

F 

T 

T 

AW* = A v(Xmin2) 

T 

F 

F 

AW* = Av(Xminl) 

T 

F 

T 

Not possible 

T 

T 

F 

AW* = Av(Xminl) 

T 

T 

T 

From Eq. 5 


If there are several constraint sets, the previous algorithms will be called 
several times for each set. The global minimum Aw will be the minimum of all of 
the solutions for each set. Additionaly, if the constraint set C x contains negative 
and positive numbers, C x will be subdivided into two different sets containing 
only negative and only positive numbers and the optimization algorithm will be 
called twice. 

Once we have an algorithm to compute the global optimal when constraints 
in x are present, all we need to do is to express different constraints of the 
problem (e.g. minimum radius, minimum flight-path-angle or maximum Aw) in 
terms of x. The next sections describe several methods to achieve this. 


5 Constraints in position magnitude r > r m i n 

In this section, we are going to describe how to impose a constraint in the 
position magnitude during the transfer: 

t G [ti,t 2 ], r(t) > r min 

We can translate the above equation into a set(s) of constraints for p. First, 
we are going to calculate the sets of p that make the eccentricy vector to be 
contained between rq and iq during the transfer and then we are going to 
compute the sets of p that satisfy r p > r m i n . The combination of these two sets 
will provide the final constraint set(s). 
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Given that the eccentricy vector is defined by (see Battin ref. ) 


e = Ar t + By 2 


where 


A = 


1 


sin 2 9 


-1 - 


r i 


T2 


— 1 cos 9 


B = 


1 


sin 2 9 


r 2 




r i 


Therefore we need to calculate the values of p that make A > 0 and B > 0, 


r 17*2 (1 — cos 9) r\r 2 (1 — cos 9) 

P — Pa — ~x~, P>Pb = w- 

r 2 — ri cos 9 r\ — r 2 cos 9 


p G [pa, +oo) n [pb, +oo) = [max(pA,Ps), +oo) = S p i — >■ e between ri, r 2 

Since p > 0, there are solutions only if min ^ j > cos 9. 

For x > 0 (short transfer) a periapse passage occurs during the transfer if 
p G S p i, for x < 0 (long transfer) a periapse passage occurs during the transfer 
if Pi S p i : 


x > 0 
x < 0 


x G 
x G 


\/ma tl(pa,Pb), +oo) = S n i 


— \/ma x(p A ,p B ),0 


= S , 


n2 


( 7 ) 

( 8 ) 


The next step will be to determine the values of p (and therefore x) that 
make r p > r pm i n . First, we will calculate e 2 = A' 2 +B 2 +2AB cos 9. After some 
algebraic manipulations we can obtain: 

2 = c 2 2 _ 2 (n + r 2 ) 2 

sin 2 9r\r^ (1 + cos 9) rir 2 P l + cos0 


or 

2 1 2 nPmirie . r,Pmin e f l r 2 

e = ~pP -2^-p + 2 — — “ 

a- cr cr r i + r 2 

Since i > 0, e,p define a hyperbola, we can rewrite the above equation to 
show this: 


1 / \2 2 Pmin 

^2 \P ~ Pmin e ) ~ e ~ ^2 


we can also proof by contradiction that 



7 * 1^2 \ 
n + r 2 ) 


( 9 ) 
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7 * 1 7*2 

Pmin e ^ ■ < 0 

n + r 2 

Therefore we can say that e and p define a vertical hyperbola centered at p = 
Pmin e , e = 0 where p is along the x-axis and e is along the y-axis (only the pos- 
itive branch is considered). The asymptotes are given by e asy = ±rl ( p — p m i ne ,) 
We can calculate the minimum eccentricy transfer by: 

de 2 = Q r,de = 2 c 2 2 (n + r 2 ) = Q 

dp ' dp sin 2 6r\r% (1 + cos 9) rir 2 


Pmin e — 


r\r 2 (ri + r 2 ) (1 - cos 9) 


If we define r pm i n := r m j n , we can calculate the constraint for p as follows: 


P 


1 + e 


> r 


pmxn 


— 1 > e 


the left-hand side of the previous equation defines a line with slope 


and 


the right-hand side of the equation is just the positive branch of the hyperbola 
defined above. The set of values of p make satisfy the constraint are the ones 
within the intersection of the line defined by — 1 and the hyperbola defined 

r pmin 

by e. Therefore, we need to first find the values of p (if any) at the intersection 
points. For that we’ll solve the following equation: 


P 

T pmin 


-i =. j 


' 2 => 2 - 2 * 


Pmin e . r%Pmin e ^1^2 
p + 2- 


d 2 r i + r 2 


1 



. ryPmirie 

r ~ + 1 d 2 

1 vmin m 


' pmin ^ J \ ' p^nm 

We have three possible solution sets: 


+ 1-2 


Pmi 


nr 2 


d 2 ri + r 2 


= 0 (10) 


p G. \prpl>prp2] 

P G [Prpl,+Oo] 

P G 0 
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Since we are solving Eq. 10 instead of 1 = e. It is necessary to 

check for extraneous solutions. We will assume that p rp i are the positive real 
solutions to Eq. 10 that also satisfy the original equation. P rp i is a valid solution 

p 

if _L£el_ - 1 

> 0. In terms of x we obtain: 

P pmin 

if two solutions exist x G \^p rp i, y/p rp 2 \ U [— y/p rp 2 , ~\fpPp l] = Sr p \ U 

if only one solutions exists x G \\Jp rp i, +oo] U [— oo, —^/p rp l] = Sr p i U Sr p 2 (12) 

if no solutions exist x G 0 (13) 

If we combine the conditions in 7 and 11, x should satisfy: 

x G Sr p 1 U Sr p2 <=>■ r p > r pmin or 


x ^ (Sr p 1 U S r p 2 ) and x ^ (S n 1 U SV12) < r pm j„ and no periapse passage 

Finally, if the initial or final terminal is inside the forbiden area, there are 
no solutions 

if n < Tmin or r 2 < Train HO Solutions 


6 Unrealistic solutions 


Not all the initial velocities in the hyperbola will generate a trajectory that 
connects the terminals. If Vi > yj 2^- and 71 > 0 the transfer trajectory will 
result on a hyperbola that goes to infinity. We need a criteria to avoid calculating 
trajectories in that region. First, we need to compute the values of p for e = 1. 
Using eq. 9 we obtain 


u _j_ 4 (nr 2 ) 3/2 _.„ 2 0 _ e _ D 

Ppar 1,2 — min e ± ^ sin - cos - = 

'min e 


4(rir 2 ) 3/2 . 9 
± 2c 2 Sm 2 Sm 


In this way the realistic set of transfers is composed of elliptical transfers 


P ^ [Pparl 5 Pp ar2 ] t t ^ G [ \/Ppar2i \/Ppri! \ ] bl [-^/Pparls \/Ppar2^] — U S e l 2 

and hyperbolic transfers with ^ < 0. We need to find a criteria to determine 
the flight-path-angle at iq. We can use the condition: 

Vl * lri (y c Xl c “t“ * Iri — 0 

We have two possible cases when x > 0: 

if i c • i r i = cos V’l > 0 —s ► x 2 < — has no solution — » 71 > 0 

v c cos ipi 
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S 7 b = 0 


if i c • i r i = cos 7/7 < 0 — > x 2 > 


v P 

v c cos t/>i 


S^b — 


V c COS 4>1 ' 



We have two possible cases when x < 0 : 


if i 0 • i r i = cos 7/7 > 0 — > x 2 > — — p has always a solution 

v c cos 7/7 


Srya — ( 0] 


if i c • i,-i = cos 7/7 < 0 — >■ x 2 < 


Vp 

V c COS 7/7 


*^7 a — 


V c COS 7/7 


•X ^ — 1 S ■-y b 

Therefore, the set of x that generates realistic trajectories is determined by: 

x £= S e i j U S e i 2 <£=> e < 1 , and 


x £ ( S e 7 U S e i 2 ) and a; <E (S ia U S 7 &) <= = => e > l,7i < 0 


We can obtain the same result using the equation that relates 71 and p from 
Battin for when 71 = 0 

c 7777 (1 — cos0) 

P = Pm 7 = 7- 

r 1 — 77 COS 0 77 — 77 cos 0 

Using, 


cost/>i 


cos ip[ 


2 2 

— r{ — c 

2 ric 


r| — rf — (rj + — 27777 cos$) 

2 riC 


—rf + 77 r 2 cos0 
77c 


where %l>[ is the interior angle between 17 and c. 
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Vp — rir 2 (1 — cos0) rir 2 (1 — cos0) 

v c cos ipi c cos ipi v ! — r 2 cos 9 

We can conclude that only tangential departure trasfers will exist if 

n 

ri — r 2 cos 9 > 0 — > cos 9 

T2 

Therefore, the set of x that generates realistic trajectories is determined by: 

x G= S e q U S e i 2 4=4 e < 1 , and 


x ( S e q U S e i 2 ) and a: G (S 7a U S 7 b) <= = => e > l, 7 i < 0 

After combining both conditions we will obtain a maximum of 4 sets (2 for 
x < 0 , and 2 for x > 0 ): S r i, -S' r2 , S r 3, SV4. 


7 Constraints on the departure and arrival flight- 
path-angles 

We can impose limits on the terminal flight-path-angles 71 and 7 2 by traslating 
those limtis into constraints on the parameter p (or x). First, we can relate 7$ 
to the parameter by 

cos 71 cos y 2 

P = PmC 7 . 777 or P = PmC 7 XV 

ri COS71 — r 2 cos(7i + 9 ) r 2 cos7 2 — r 1 cos(7 2 — 0) 

If 7mm A 7i *7 'Imax: 

dp 

P e b (7m«x) ,P (7mm)] for terminal 1, since — < 0 

« 7 i 

dp 

P G [p ( 7 min) ,P (7mm)] for terminal 2, since — - > 0 

« 7 2 

8 Constraints on the arrival and departure An 

We can impose constraints in the magnitude of the maneuvers A by solving 
a quartic polynomial in x. First, by using the definitions, 

A Vi = || Vj - Vqoll < A Vi, min 4=4- (v» - V i; o) • (Vj - Vq 0 ) < A v\min 

This will lead to the solution of the following quatic polynomial for terminal 

1: 
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Table 5: Characterization of the solutions to the maximum A Vi contraint prob- 
lem. 


real solutions to Eqs. 14 and 15 : 2 q < x 2 < £3 < 24 

x > 0 x < 0 


0 solutions — > C xl = 0 

2 solutions — ► C x2 = [xi, x 2 ) 

2 solutions — > C x 1 = [cci,a; 2 ] 

0 solutions — > C x2 = 0 

2 solutions — ► C^i = [^3,2:4] 

2 solutions — > C x2 = [aq, x 2 \ 

4 solutions ->• C xl = [a? 3 , ar 4 ], C x2 = [£3,24] 

0 solutions 

0 solutions 

4 solutions -> C x 1 = [£3,24], C x2 = [24, 24] 

0 solutions— > C x \ = 0 

0 solutions — > C x2 = 0 


v c x 4 -v c (i c • V10) x 3 +[ 2 v c v p (i c • i r i) + v\ 0 - Avf min ] x 2 - 2 v p (i r i • vi, 0 ) x+v 2 p = 0 

(14) 

and for terminal 2: 

v c x 4 —v c (i c • v 2>0 ) x 3 +[~ 2 v c v p (i c • i r2 ) + vl 0 - Av 2 min ] x 2 + 2 v p (i r2 • v 2j0 ) x+v 2 = 0 

(15) 

Given the signs of the polynomials in Eqs. 14 and 15, we can establish the 
nature of the solutions (see Table 5) and from them we can create the constraint 
sets. 
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